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Computer modeling of multicellular systems has been a valuable tool for interpreting and guid- 
ing in vitro experiments relevant to embryonic morphogenesis, tumor growth, angiogenesis and, 
lately, structure formation following the printing of cell aggregates as bioink particles. Computer 
simulations based on Metropolis Monte Carlo (MMC) algorithms were successful in explaining and 
predicting the resulting stationary structures (corresponding to the lowest adhesion energy state). 
Here we present two alternatives to the MMC approach for modeling cellular motion and self- 
assembly: (1) a kinetic Monte Carlo (KMC), and (2) a cellular particle dynamics (CPD) method. 
Unlike MMC, both KMC and CPD methods are capable of simulating the dynamics of the cellular 
system in real time. In the KMC approach a transition rate is associated with possible rearrange- 
ments of the cellular system, and the corresponding time evolution is expressed in terms of these 
rates. In the CPD approach cells are modeled as interacting cellular particles (CPs) and the time 
evolution of the multicellular system is determined by integrating the equations of motion of all 
CPs. The KMC and CPD methods are tested and compared by simulating two experimentally well 
known phenomena: (1) fusion of two spherical aggregates of living cells, and (2) cell-sorting within 
an aggregate formed by two types of cells with different adhesivities 

PACS numbers: 87.17.Aa, 87.17.Rt, 87.85.G-, 87.85.Lf 
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I. INTRODUCTION 

Understanding how living cells form tissues and organs 
is a fundamental problem of developmental biology pj [2j , 
and is also important for the rapidly expanding field of 
tissue engineering that aims at building functional tis- 
sue substitutes in vitro [3]. Tissue engineered structures 
may be used for drug testing and to restore or replace 
damaged tissues and organs [4] . An emerging tissue engi- 
neering technique is bioprinting [5HTT] via the automated 
layer-by- layer deposition of multicellular aggregates (the 
bioink). Subsequent postprinting fusion of the contigu- 
ous aggregates gives rise to the desired tissue construct. 
Predicting the result of post-printing tissue formation is 
a task for theoretical modeling. 

A guiding principle for most models of cell rearrange- 
ment in cell aggregates is the differential adhesion hy- 
pothesis (DAH) proposed by Steinberg [12l [13]. DAH 
states that structure formation in multicellular systems 
occurs due to (i) differences in cell-to-cell adhesion of 
different types of cells and (ii) cell motility. Cells seek 
positions with the largest number of strong bonds. For 
example, in a random mixture of two cell types of differ- 
ent cohesivities the more cohesive cell population sorts 
out and occupies the central region surrounded by the 
less cohesive population. 
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By incorporating DAH, Metropolis Monte Carlo 
(MMC) simulations correctly predict the formation of 
multicellular structures of minimum energy of adhesion, 
and identify long-lived, metastable configurations [14J . 
However, MMC cannot predict the actual time evolution 
of multicellular systems. 

Insight into time evolution emerged from experiments 
designed to verify DAH, which revealed that embryonic 
tissues behave analogously to highly viscous liquids. The 
concept of tissue liquidity motivated a quantitative de- 
scription of embryonic tissues in terms of parameters of 
continuum hydrodynamics, such as surface tension and 
viscosity. For example, the (apparent) surface tension 
7 was measured for several tissue types using a paral- 
lel plate compression apparatus [15 and the values were 
found to be consistent with the sorting behavior of these 
tissues [16 . 

Here, by using the concept of tissue liquidity and DAH, 
we formulate two alternatives to the MMC approach for 
modeling cellular motion and self-assembly: (1) a kinetic 
Monte Carlo (KMC), and (2) a cellular particle dynam- 
ics (CPD) method. Both these methods are capable 
of describing and predicting the real time evolution of 
the shape of multicellular systems and tissue constructs 
in certain morphogenetic processes and post-bioprinting 
structure formation. In the KMC method the configu- 
ration of the multicellular system is propagated in time 
through a standard rejection-free kinetic Monte Carlo al- 
gorithm. This approach should provide a more accurate 
description of the time evolution of a multicellular system 



than other grid based methods, such as, the MMC model 
[61 [H] or the widely used cellular Potts model (CPM). 
The latter uses a modified Metropolis MC algorithm to 
update the configuration of the simulated system and 
postulates that time is proportional to the number of 
MC steps, which in general is not the case [17]. In the 
CPD method [17 individual cells are modeled as an en- 
semble of cellular particles (CPs) that interact via short 
range contact interactions, characterized by an attractive 
(adhesive interaction) and a repulsive (excluded volume 
interaction) component. CPs in a cell are held together 
by an additional confining potential that mimics the role 
of the cell membrane. The time evolution of the spatial 
conformation of the multicellular system is determined 
directly by recording the trajectories of all CPs by inte- 
grating their equations of motion. What sets apart CPD 
from the other similar off- grid particle methods, such as 
Newman's subcellular element method (SEM) [TSHST] . is 
the employed force field (especially the confining poten- 
tial) and its parametrization that makes the system be- 
have as a complex viscous liquid. In particular, the CPD 
model parameters are determined such that the shape of 
two fusing spherical aggregates in the CPD simulation 
match as closely as possible the one observed experimen- 
tally, i.e., two attached spherical caps (see Fig. [l]) [22] . 
To test and compare the KMC and CPD methods, we 
apply them to simulate the fusion of two spherical aggre- 
gates and the evolution of cell sorting within an aggre- 
gate, two morphogenetic processes driving postprinting 
structure formation. For the theoretical description of 
the fusion of two identical spherical aggregates we use a 
simple continuum model introduced by Frenkel [23 and 
further developed by others working in the field of rhe- 
ology [24| [25] . It is this theoretical continuum model 
that provides the link between the time scales of simula- 
tions and the time scales of experiments. Once this link 
is established, the KMC and CPD simulations are used 
to quantitatively predict the time evolution of complex 
postprinted structures whose description using a contin- 
uum hydrodynamics approach is impractical. 

The remainder of the paper is organized as follows. 
Secti on [ll| describes the KMC (Sec. |IIA| ) and the CPD 
(Sec. II B) methods, as well as the theoretical asp ects o f 
the continuum approach of aggregate fusion (Sec. II C). 
Section [1111 contains the results and discussion of our 
KMC and CPD simulations, i.e., f usion o f identical spher- 
ical multicellular aggregates (Sec. Ill A) and cell sorting 



approaches equilibrium, or is in a metastable state, the 
Metropolis algorithm rejects most trial moves because 
the acceptance probability is small. A main feature of 
the KMC algorithm is that it is "rejection- free" . In each 
step, one calculates the transition rates for all possible 
changes compatible with the current configuration, and 
then chooses a new configuration with a probability pro- 
portional to the rate of the corresponding transition. 

We designed and implemented a KMC algorithm to 
simulate the time evolution of a lattice model of multicel- 
lular systems. Aggregates of cells in cell culture medium 
are represented on a 3D hexagonal close-packed lattice 
by associating each site to either a cell or to a similar 
sized volume element of medium. Thus, the lattice spac- 
ing is equal to one cell diameter. We assume that each 
cell interacts with its 12 nearest neighbors (1st and 2nd 
neighbors considered to be nearest neighbors) located at 
a distance of one lattice spacing from the given cell. Inter- 
actions are expressed in terms of works of cohesion and 
adhesion [27l [28] , defined as the work needed to break up 
the contact between two neighbors of respectively simi- 
lar or differing types of cells. For example, in case of 
a multicellular aggregate composed of a single cell type, 
the work needed to extract a cell from the aggregate (i.e. 
model tissue) is the work of cohesion, Ccc^ multiplied by 
the number of the cell's nearest neighbor. The interac- 
tion between cells and the cell culture medium is set to 
zero. The movement of cells is described by assigning 
rates to swapping cells with adjacent cells of different 
type and/or with medium elements. These elementary 
moves occur with rates given by 
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-E^/Ej 



(1) 



(Sec. IIIB). Conclusions are presented in Sec. IV 



where the factor wq is the frequency of attempts to cross 
the energy barrier of height ^^b, and Et is the energy 
of biological fiuctuations [29 , the analog of the energy 
of thermal fiuctuations, ksT (ks is Boltzmann's con- 
stant and T is the absolute temperature). It has been 
argued that Et is a characteristic measure of cell motil- 
ity: the higher is Et in comparison to the energies of 
cohesion/adhesion, the higher is the motility of the cell 
[291. 

Due to the complexity of the cytoskeletal machinery 
responsible for cell movement, there is no unique way 
to assign a barrier height to the swapping of two cells. 
Any reasonable choice, however, needs to be consistent 
with the following set of experimental observations on 
cell movement in 3D: 



II. COMPUTER AND THEORETICAL 
MODELING 

A. Kinetic Monte Carlo for Multicellular Systems 

The Kinetic Monte Carlo method (KMC) was pro- 
posed as an alternative to the MMC method for simu- 
lating the evolution of Ising models ^26^ . When a system 



(1) Relocation of cells in embryonic tissues and in some 
engineered tissues (such as cell aggregates) occurs 
according to DAH [l3", 30 : cells take advantage of 
their motility to establish the maximum number of 
strong bonds with their neighbors. 

(2) Anchorage-dependent cells do not spontaneously 
dissociate from the cell aggregate they are part of 



(3) The speed of cell movement in 3D matrices has 
a particular dependence on the strength of cell- 
matrix adhesion: cell movement is fastest at an op- 
timal strength of binding. Too weak or too strong 
binding hampers cell movement [31^ .32^ . 

Consider a binary particle model for a multicellular 
system formed by two cell types, t = 1,2 (for a multi- 
cellular aggregate composed of one cell type, t = 1, sur- 
rounded by tissue culture medium, t = 2 represents the 
medium particle). The configurational energy (or total 
interaction energy), E is expressed as [6] 



E = 7i2A^i2 + const, 



(2) 



"m" for which Km-i < uKm < Km] (S5) Generate an- 
other uniform random number u' between and 1, and 
increment the time variable (i.e., t ^ t ^ At) by the 
non-uniform time step 



At = -K^^ log(w'); 



(5) 



(S6) Update all rates kn that may have changed due to 
the previous transition "tti" ; (S7) Return to step S2 and 
repeat the process until the time variable reaches the 
desired target value. 






where 712 = (en + e22)/2 - 612, with en and 622 being 
the energies of cohesion respectively for cell type 1 and 

2, and 612 is the energy of adhesion. A^i2 

Si=i ^ii is total number of nearest neighbor pairs of 
different cell types cells, n^2 (^ii) the number of nearest 
neighbors of cell i of type 1(2), which are of type 2(1) and 
NI {N^^) the total number of cells of type 1(2), which 
have at least one (nearest) neighbor of type 2(1). (As the 
const is irrelevant for the evolution of the system, we set 
it to zero [&|.) 

Consider two nearest neighbor cells, i and j of different 
types (without loss of generality we can set z = 1 and j = 
2). The system evolves in time towards configurations of 
decreasing energy E, i.e. for 712 > (712 < 0) N12 
decreases (increases). For 712 > and 712 < cells 
respectively phase separate (cell sorting) and mix (cell 
mixing). Elementary KMC moves consist of swapping 
two neighbors of different types (swapping cells of same 
type does not change the energy). The contribution of 
two such cells, i and j to E is 



^. 



ij - -{ni2^nji), 



(3) 






and E = ^i2i^j=i^ij' Furthermore, the larger is Eij 
the more likely is the KMC move to swap cells i and j. 
Thus it is reasonable to define the energy barrier E^^ in 
Eq. ([1]), for a transition involving the swapping of two 
cells i and j, as 



0<El^ = E^^^ 



Ej, 



(4) 



where E^^^ is the maximum possible value of Eij . For 
712 > 0, E^^^ is obtained when the number of neighbors 
of differing type surrounding cells i and j is maximal. 

Now we can formulate the steps of our KMC algo- 
rithm for simulating the time evolution of multicellu- 
lar systems: (SI) Set t = 0; (S2) Find all interfacial 
cells (i.e., cells in contact with cell culture medium or 
with cells of different type) and compute the rates /c^, 
1 < m < M, corresponding to all possible M transitions 
involving these cells; (S3) Calculate the cumulative rates: 
Km = Xl!r=i ^n7 1 ^ ^^ ^ M; (S4) Generate a uniform 
random number u between and 1 and carry out event 



B. Cellular Particle Dynamics Method for 
Multicellular Systems 



The cellular particle dynamics (CPD) method is an 
off-lattice, particle-based computer simulation method 
that can describe and predict the time evolution of 3D 
multicellular systems during shapechanging biomechan- 
ical transformations [17]. Within the CPD formalism 
cells, regarded as continuous objects with self-adaptive 
shape, are coarse-grained into a finite number, Nqp^ of 
equal volume elements. Each volume element is repre- 
sented by a point-like cellular particle (CP) situated at 
its center of mass. CPs interact via short-range contact 
interactions, characterized by an attractive (adhesive in- 
teraction) and a repulsive (excluded volume interaction) 
component. In addition, CPs within a given cell are sub- 
ject to a confining potential that assures the integrity of 
the cell. The time evolution of the spatial conformation 
of the multicellular system is determined directly by cal- 
culating the trajectories of all CPs (and, therefore, cells) 
through integration of their overdamped Langevine equa- 
tions of motion. This minimalist model, when properly 
parametrized, has the features of a complex viscous liq- 
uid and it is suitable for describing the time evolution of 
multicellular aggregates and soft-tissue constructs. 

For the n}^ CP in cell a, the equation of motion is 



/ir,.(t) = -Vc.^/7 + /«Jt), (6) 

where r^^ {t) is the position vector, U is the potential en- 
ergy function describing the interaction of the CPs , fi is 
the friction coefficient, /^^(t) is a random force, and the 
dot denotes time derivative. The random force is mod- 
eled as a Gaussian white noise with zero mean and vari- 
ance {fi{t)fj{0)) = 2Dii^5{t)5ij, where D is the sort-time 
self diffusion coefficient of the CPs. The CPD parame- 
ters D and ji are related to the previously introduced 
biological fluctuation energy Et by the Einstein relation 
Dfi = Et. The CP potential energy U has an intra- 
cellular and an inter-cellular component corresponding 
to CPs belonging respectively to the same cell and to 



different cells, i.e., 






■lEE^' 



■rpj), (7) 



where an (Pm) labels the cellular particle n {m) in cell a 
(P). We model the short-range contact inter- and intra- 
cellular interactions between CPs through 

W'^'^^'^ir) = Vlj (r; e^^'^^ a^^'^^) , (8a) 



W''''^''{r) = Vlj (r; e^^'^^ a^^'^'^) + - (r - 0^ ©(^ - 0, 



where 



(8b) 



VLjir;e,a) = 4e 



(;)"-©' 



is the standard Lennard- Jones potential, and 0(r) is the 
Heaviside step function. Note that instead of Vlj one 
could use any other potential that has a repusive core 
and a short range attractive part. For example, in SEM 
[18] a Morse potential is used to describe the interaction 
between two subcellular elements. However, the impor- 
tant addition in CPD is the second quadratic term in 
jjmtra ^j^^^^^ for r > ^, represents an elastic confining po- 
tential used to maintain the integrity of the cell. This 
term, characterized by the elastic constant /c, guarantees 
that the CPs within a cell remain confined inside the 
boundary of the cell. The time evolution of the multi- 
cellular system within the CPD approach is determined 
by numerically integrating the equations of motion (|6| 
for all CPs. We have accomplished this by implementing 
the intra- and inter-cellular interaction forces, Eqs. ([t])- 
(Is]), and a Langevin dynamics integrator in the freely 
available massively parallel molecular dynamics package 
LAMMPS [33]. 



1. CPD units and parameters 

For a multicellular system of a single cell type there are 
nine CPD model parameters that need to be determined: 
Ncp (the number of CPs per ceh), D, /i, a^^^^^, e^^^^^, 
/c, ^, cr*^^^^ and e*^^^^. The choice of Ncp is determined 
by the degree of detail we want to describe individual 
cells. Since we are interested in the time evolution of 
the shape of an aggregate formed by a large number of 
cells and not in the detailed description of the surface 
dynamics of individual cells, in the present work we make 
the reasonable choice of Nqp = 10. Because a in (8c) 



determines the size of the interacting CPs, we can set 



a = a 



winter 



The length ^ in (8b) represents the 



size (diameter) of a cell, which comprises Ncp tightly 

1/3 

packed CPs of size a. Thus, one can estimate ^ ~ ^^cp • 



Next, we define convenient CPD (or computer) length, 
time and energy units according to 



^0 = cr, to = 



D 



Eo = Et = i^D, 



(9) 



In these units all CPD parameters all pure numbers, and 
in particular a = D = /j. = 1. We set the confining poten- 
tial parameters as ^ = 2.5 (~ 10^/^) and k = 5. A larger 
(smaller) value for k makes the cell more rigid (soft) when 
subjected to deformations. The chosen value, for which 
ka'^/2 = 2.5Et, is suitable when cells in the aggregate 
are exposed only to adhesion and surface tension forces. 
By choosing e*^^^'^ = 1 (i.e., same as the biological fiuc- 
tuation energy Et) the dynamics of the CPs inside a cell 
will have sufficient randomness to produce cell surface 
ffuctuations that play an important role in cell motility 

124. 

Thus, out of the nine CPD parameters we are left with 
only one, e'^^^^^ ^ that needs to be determined such that 
the time evolution of the shape of the multicellular sys- 
tem follows as closely as possible the corresponding ex- 
perimental one. For this purpose we focus on the fusion 
of two identical spherical aggregates, as described in the 
next section. Based on extensive CPD simulations we 
have found that the best agreement with experiment is 
obtained for 1 < e*^^^^ < 2, when the system behaves as 
a viscous liquid. By increasing e*^*^^ above 2 the fusing 
cellular aggregates show sign of solidification and their 
behavior deviate significantly from experiment. The re- 
sults reported in this paper are for e*^^^^ = 1. However, 
these are similar to the ones obtained for any e*^*^^ < 2. 

Furthermore, in all our CPD simulations we have used 
an integration time step At = 10~^, and used a cutoff 
radius Re = 2.5 for [/^^^^^(r). 



C. Continuum Description of the Fusion of Tavo 
Spherical Cell Aggregates 

The fusion of two contiguous cell aggregates is driven 
by surface tension, 7, and resisted by viscosity, 77. It is 
an experimental fact that during the fusion of identical 
spherical soft tissue aggregates the shape of the system is 
that of two touching spherical caps (see Fig.fTl) [22 . This 
observation suggests that soft tissues behave like complex 
viscous liquids whose description requires an a priori un- 
known hydrodynamic constitutive model. However, the 
simplicity of the geometry allows us to describe analyt- 
ically the dynamics of the considered fusion process by 
employing conservation laws as proposed by Frenkel [23] 
and Eshelby [35] for the coalescence (sintering) of highly 
viscous molten drops. 

The fusing aggregates are modeled as two spherical 
caps of radius R{0) with circular contact ( "'neck'''' ) region 
of radius r{9) = R{9)^\ii9 (see Fig. [l|\). Volume conser- 
vation requires 

R{e) = 2^/^(1 + cos6>)-^/^(2 - cos6>)-^/^i?o , (10) 




Note that according to Eqs. (|T2|) and (14) the dynamics 



FIG. 1. (color online). (A) The shape of two fusing spherical 
cell aggregates can be quantified by the angle 0{t) and the 
radius R(t). (B) Throughout the fusion of two indentical 
spherical cushion tissue aggregates [22 the system has the 
shape of two connected spherical caps. The numbers indicate 
the time (in minutes) elapsed from the start of the fusion. 



with Rq = R{0). Thus, the time evolution of the fu- 
sion process is parametrized by a single angle = 6{t)^ 
defined in Fig. [TJA, that changes from 0{0) = to 
^(oo) = 7r/2. The rate of the decrease in surface en- 
ergy is Ws = jdS/dt^ where the free surface area S = 
S{e) = 47ri?2(l9)(l + cos l9). The equation of motion for 
6{t) can be derived by equating Ws with the rate of the 
energy dissipated by the viscous flow VK^ ^ —ATrRlrja^ 
[23[ l25]. Assuming biaxial stretching flow, 



dvx 
dx 



1 d 

'R(0)Jt 



[R{0) cos 0]. 



(11) 



Inserting Eq. ( 11 ) into the energy balance equation Ws 



Wrj leads to 



dO _ 1 sin (9 cos (9(2 - cos 19)^/^ 

dt ~ r 2^/3(1 -cos6>)(l + cos 6>)V3 ~ 2^ R{0) 

where the characteristic fusion time 



1 RocotO 



(12) 



r = r]Ro/j . 



(13) 



Equation (12) can be solved numerically for = 0(t). 



However, one can derive a simple and accurate analytical 
approximation for 9{t) by setting R{9) ^ Rq in Eq. ( p!2| ). 
Indeed, throughout the fusion process 1 < R{0)/Rq < 
2^/^ ^ 1.26 holds. With this approximation, Eq. ([l2| 
can be easily integrated with the result 



of the fusion process, described by 6>(t)as a function of 
t/r, is independent of the size (i.e., Rq) of the fusing 
spheres. Rq appears only in the characteristic fusion time 
r, Eq. ([13|. 



Finally, using Eq. (14), the square of the radius of the 



circular neck region of the fusing spherical caps can be 
expressed as 



r 
Rq 



A(t)[l-exp(-t/r)], 



(15a) 



with 



A{t) = 2^/3 (l + e-*/2- j '^' (2 - e-*/2-) '^' . (15b) 

For short times, t ^ r, Eqs. dTsl yield the familiar 
linear-in-time expression, (r/Ro) ~ t/r, obtained by 
Frenkel [23 and Eshelby [35 . While this formula has 
been applied previously to estimate the capillary veloc- 
ity Vc = j/t] = Rq/t of soft tissues [22", "ST, we are not 
aware of any previous study that followed the time evolu- 
tion of the shape and of [r(t)/i?o]^ throughout the fusion 
process of two spherical tissue aggregates. First, we have 
determined the dimensionless CPD parameters (i.e., ex- 
pressed in CPD units; see Sec. IIB 1 ) such that the shape 



of the fusing aggregates during CPD simulation resemble 
as close as possible to spherical caps. Second, we have 
determined the characteristic fusion time r by fitting the 
data for [r{t)/ Rq]^ ^ obtained respectively from experi- 
ment, CPD and KMC simulations, to Eqs. (15). Finally, 
the CPD time unit can be calculated as to = Texp/TcPD- 
Once to is known, one can predict through CPD simu- 
lation the time evolution of an arbitrary 3D tissue con- 
struct built from the same type of cells for which the time 
calibration was performed through the above method 
(i.e., fusion of spherical aggregates). Clearly a similar 
calibration strategy can be used for the KMC method. 



III. RESULTS AND DISCUSSION 

To test and compare the KMC and CPD methods 
described in Sec. [llj we have applied them to simulate 
two important morphogenetic processes: (A) tissue fu- 
sion (the fusion of two identical spherical multicellular 
aggregates), and (B) cell-sorting (within a spherical mul- 
ticellular aggregate formed by two types of cells with dif- 
ferent adhesivities). 



A. Fusion of Two Spherical Cell Aggregates 



As descirbed in Sec. |IIC[ the fusion of two identical 
spherical aggregates can quantitatively be characterized 
by the time dependence of the radius, r(t) , of their circu- 



cos^ = exp(— t/2r) 



lar contact region. According to Eqs. (15), r(t) obtained 



(14) from experiment and from KMC and CPD simulations. 



can be used to determine the characteristic fusion time 



r, Eq. (13). Thus, for a given ceh type, by comparing 
the experimental r with that obtained from computer 
simulations one can calibrate the time scale of the cor- 
responding computer model. Once such a calibration is 
done, one can make quantitative in silico predictions of 
the time evolution of various multicellular processes that 
involve the same cell type [17]. 

In this section we present KMC and CPD simulation 
results for the fusion of two identical spherical aggregates. 
We show that in both cases the computed (r/i?o)^ vs t/r 
dependence can be reasonably well fitted by Eqs. (15). 
Then, using experimental results for aggregate fusion 
[22] . the calibration of the KMC and CPD simulation 
time scales is exemplified for the case of cardiac cushion 
tissue (CT). Finally, KMC and CPD simulations are used 
to predict the formation of a toroidal structure by cell ag- 
gregate fusion, an important structure in the engineering 
of tubular tissue constructs 8 . 



1. KMC simulations 

The initial radius of the two identical fusing aggregates 
used in our KMC simulation was Rq = 10 cell diame- 
ters. Each aggregate contained 5,927 cells, with a cell-cell 
work of cohesion ecc = 0.9. The medium- medium (cell- 
medium) work of cohesion (adhesion), emm (ecm), was 
considered to be negligibly small. A total of 10 KMC 
simulations of the same fusion process were carried out, 
each time using a different seed of the random number 
generator. Each simulation was run for 10^ KMC time 
steps. 

Representative snapshots during the KMC fusion sim- 
ulation are shown in Fig. ^ The corresponding (r/i?o)^ 
vs t/r dependence is shown in Fig. [s] (dashed curve). 
Apart from the beginning of the fusion process (i.e., 
t < r) the KMC result appears to match rather well both 
the theoretical prediction (thick-solid curve), Eqs. ( [T5| ), 
and the experimental results corresponding to the fusion 
of CT aggregates (open-circle) [22 . 

In the KMC time unit to = Wq , the fusion time [ob- 
tained by fitting the KMC simulation results to the the- 
oretical formula Eqs. (15)] was tq = 1.1 x 10^. Since 
the experimental characteristic fusion time for CT ag- 
gregates Texp ^ 5h [22 , it follows that the KMC time 
unit (for CT aggregates used in [22 ) has the calibrated 



KMC 



CPD 



value to 



Texp/ro = 1.6 X 10 ^s. 



2. CPD simulations 



Each of the two spherical aggregates used in the CPD 
simulation of aggregate fusion contained 2000 cells. The 
used CPD parameters and integration timestep, in CPD 
units, are described in Sec. |IIB 1[ The equilibrated ag- 




FIG. 2. Time evolution of the fusing aggregates in the KMC 
(left) and CPD (right) simulations. The snapshots were taken 
at: (a) t = 0, (b) t = O.lQr, (c) t = 2.8t, and (d) t = 5.5t. 
The solid-line contours represent the theoretical shapes of the 
fusing aggregates determined by Eqs. (15). 



Representative snapshots during the fusion process are 
shown, and compared with the corresponding KMC sim- 
ulation results, in Fig. [2] While in both KMC and CPD 
simulations the profiles of the fusing aggregates for in- 
termediate stages of the fusion process (Fig. [2b-c) agree 
quite well, these show noticeable differences with respect 
to the theoretical prediction, Eqs. ( [T5| ), shown as solid- 
line contours in Fig. [2] 

The (r/i?o)^ vs t/r dependence in the CPD simula- 
tion is also shown in Fig. [3] The CPD and KMC simula- 
tion results are similar. Apart from short times {t < r) 
they agree quite well with both the theoretical prediction, 
Eqs. (15), and the experimental results for CT [22] , 



gregates were placed within a distance of one a before 
starting the fusion simulation. 



The characteristinc CPD fusion time is determined to 
be r ~ 540^0. By equating this with Texp ^ 5h, one finds 
that the CPD time unit calibrated for CT aggregates is 
to ^ 0.6 min. The CPD simulations were preformed on 32 
CPUs of a dual core 2.8GHz Intel Xeon EM64T cluster 
with a performance of around 5 million timesteps/day 
(which is equivalent to 500to and slightly less than Ir). 
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FIG. 3. (color online). Comparison of (r/Ro)'^ vs t/r for 
the fusion of two spherical aggregates obtained from KMC 
simulations (dashed line), CPD simulations (thin solid line), 
continuum theory (thick solid line) and experiment (circles) 
using cardiac cushion tissue (CT) aggregates ^22j. 



3. Toroidal structure formation 

Once the KMC and CPD time scales have been cali- 
brated from the fusion of two spherical CT aggregates, 
one can employ KMC and CPD simulations to describe 
and predict the time evolution of more complex CT struc- 
tures, which are not tractable analytically. To exemplify 
this point, here we consider the formation of a toroidal 
structure as a result of the fusion of 10 identical CT 
spherical aggregates initially arranged in a circular con- 
figuration as shown in Fig. |4^. The corresponding KMC 
and CPD simulations were carried out using the same 
model parameters as in the fusion of two aggregates de- 
scribed above. In both KMC and CPD simulation the 
fusion process into a toroidal ring appeared to be com- 
pleted in At ^ 2.5r ^ 12.5 h, as shown in Fig. (4)3. This 
prediction should be easily testable experimentally. 

While it seems that both KMC and CPD methods are 
capable of providing a fairly good description of the shape 
evolution of a multicellular system during its biomechan- 
ical relaxation process, the actual cellular dynamics in 
the two methods is quite different. Indeed, unlike in 
CPD simulations, in KMC simulations the motion of in- 
dividual cells is unrealistically fast. This point is man- 
ifest in Fig. [4j By the time the toroidal ring structure 
is formed, in the KMC simulation, cells from adjacent 
aggregates (colored differently) appear to be completely 
mixed. This is clearly not the case in the CPD simu- 
lations, where, similarly to existing experimental results 
[HI [22] , there is little mixing between the cells of the fused 
adjacent aggregates. 

To further emphasize this point, we have quantified 
the degree of cellular mixing during the fusion, along the 
X-axis, of two identical spherical aggregates [labeled as L 
(left) and R (right)], with initial radius Rq (see Fig. fl]). 




FIG. 4. KMC (left) and CPD (right) simulations of toroidal 
structure formation through the fusion of 10 cell aggregates. 
Top view of the fusing aggregates at (a) the beginning (t = 0), 
and (b) the completion of fusion, (c) Cross-section through 
the median plane of the fused toroidal structure shown in 
(b). Otherwise identical cells, initially located in adjacent 
aggregates are colored differently to emphasize the degree of 
mixing during fusion. 



by calculating the time dependent mixing parameter 
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[A7Vm(t)]2 • 



(16) 



Here AN^{t) [A^{t)] is the number of CPs situated ini- 
tially (at t = 0) in the L (R) aggregate and having, at 
time t, the x{t) coordinate in the interval {—2Ro -\- {m — 
l)Ax, -2Ro + mAx}, 1 < m < M, with M a prop- 
erly chosen, sufficiently large integer. Ax = ARq/M, and 
ANm{t) = AN^{t) + AN^{t). Clearly, dmix can take 
values between (completely unmixed system) and 1 
(uniformly mixed system). 

The time evolution of dmix{t) is shown in Fig. [s] In 
the KMC simulation cellular mixing is almost complete 
{dmix = 1) after the characteristic fusion time r, i.e., sig- 
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FIG. 5. Time evolution of the mixing parameter dmix calcu- 
lated for the fusion of two cellular aggregates from the CPD 
(solid line) and KMC (dashed line) simulations. 



nificantly sooner than the completion of the fusion pro- 
cess (~ ^t\ By contrast, in the CPD simulation even at 
the end of the fusion dmix ^ 0-2 <C 1. Based on these 
results one may conclude that: (i) the cellular dynam- 
ics that drives aggregate fusion in the KMC simulations 
is unrealistic (i.e., the system is too liquid-like), and (ii) 
the CPD model provides a more realistic and atractive 
approach to describe biomechanical relaxation processes 
of multicellular systems. 



B. Cell Sorting in Tavo Component Aggregates 



TABLE I. Values of the model parameters (energies expressed 
in units of Et^ used in the KMC and CPD simulations shown 
in Fig. [6] 
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cal aggregate into two attached homogenous spheroidal 
caps (each containing either a or 6 cells) as shown in 
Fig. pp. Thus, the degree of cell sorting is enhanced (re- 
duced) for small (large) values of the adhesion energy 
6^6, compared to the corresponding cohesion energies e^a 
and 655. Note that in terms of the interfacial tension jab 
(defined below Eq. ([2|, for "1" = a and "2" = 6), case 
CI corresponds to 7^5 > and Cab > ^aai while case C2 
corresponds to 7^5 < 0. The inequalities defining case 
C3 also imply 7^5 > 0. Thus, in a multicellular aggre- 
gate with two types of cells, in order to have cell sorting 
(segregation) the corresponding interfacial tension must 
be positive (i.e., 7^5 > 0). The larger this parameter the 
more efficient and complete the sorting. 

The results of our KMC and CPD simulations, pre- 
sented next, appear to be in good agreement with in 
vitro experimental findings for these three cases [28,. 



When two populations of cells of different adhesivi- 
ties are randomly mixed within a multicellular aggregate, 
they sort such that the more adhesive cells occupy the in- 
ternal region while being surrounded by the less adhesive 
cells. Cell sorting has been extensively studied both in 
vitro fT^[36H3H] and in silico [39-41 . 

According to DAH, the outcome of cell sorting in 
a two-component multicellular aggregate (composed of 
two types of cells, labeled 'a' and '6') depends on 
the relative magnitude of the corresponding works 
of cohesion/adhesion needed to separate cells of the 
same/different types (i.e., e^a, ^55, and 6^5), respectively 
[28] . Here we employ both KMC and CPD simulations 
(described in Sees. II A and II B) to investigate cell sort- 



ing in a spherical aggregate of two cell types a and 6, with 
^aa < ^bb' We cousidcr three cases, referred to as CI, C2 
and C3, that lead to qualitatively different experimental 
outcomes [28 . CI: For intermediate adhesion between a 
and b cells, i.e., Caa < ^ab < (eaa + e55)/2, the less cohesive 
a cells engulf the more cohesive b cells, thus leading to 
the complete segregation (see Fig. |6b). C2: For strong 
a-b adhesion, i.e., (e^a + ^bb)/'^ < ^ab-, there is limited 
sorting and the spherical aggregate remains more or less 
homogeneously mixed (see Fig. JGJ^). C3: For weak a- 
b adhesion, i.e., 6^5 < ^aa < ^bbi the two types of cells 
completely separate by transforming the initial spheri- 



1. KMC simulations 

We have performed three KMC simulations of cell sort- 
ing starting with a spherical aggregate composed of a 
random mixture oi Na = 3, 589 less cohesive cells of type 
a and A^5 = 2, 362 more cohesive cells of type b (i.e., with 
^aa < ^bb)' Thus, the Spherical aggregate had a total 
of A^ = 5, 951 cells, and a radius of about 10 cell diam- 
eters. The values of the model parameters used in the 
three KMC simulations, corresponding to cases CI, C2 
and C3 described above, are listed in Table [l| Each KMC 
simulation was performed up to 10^ (non-uniform) time 
steps, given by Eq. ([5|, leading to the final configurations 
shown in Fig. ^p-d. 

To quantify the degree of cell sorting as a function 
of time during the KMC simulations, we used a sorting 
parameter s defined as [43j 



N ^ Ni 

i—1 



(17) 



where N is the total number of cells in the system, and for 
a given cell z, Ni (Nt-) is the number of nearest neighbor 
cells regardless of their type (of the same type ti as the 



KMC 



CPD 



the interface NtjNi ~ 1/2, according to Eq. (17) 




FIG. 6. 3D snapshots from KMC (left) and CPD (right) sim- 
ulations of the cell sorting in an initially spherical aggregate 
composed of two randomly mixed cell types (black and light 
grey). The snapshots represent the (a) initial, and (b-d) final 
configurations of the simulated system. The latter correspond 
to (b) intermediate (case CI), (c) strong (case C2), and (d) 
weak (case C3) cell adhesion energy, as explained in the text. 
For better visualization of cell mixing/sorting in (a)-(c) only 
half of the spherical aggregate is shown. (Images rendered 
with VMD Ha). 



cell i). The sum in Eq. (17) runs over all cells in the sys- 
tem. Clearly, < 5 < 1, and the larger s the more com- 
plete the sorting. Note that even for completely sorted 
multicellular systems, built from two (or more) different 
cell types, the presence of the interface (s) between the 
segregated regions renders the maximum possible value, 
Smaxi of the sorting parameter Smax < 1- For example, 
in the above case CI, when at the end of sorting Na cells 
of type a completely engulf Nij cells of type 6, one can es- 
timate Sjnax as follows. For simplicity, assume that both 
cell types have spherical shape with the same diameter 
d. Let AA^ be the number of cells (of either type a or 
b) situated at the spherical interface, of mean radius R^ 
and width Ai?, between the two segregated regions (see 
Fig. pp), and N = Na + Ni^. Since for a cell i situated at 



TV 



- X AA/- + 1 X (TV - ATV) 



IAN 



uni- 

^6' 



Furthermore, assuming that cells are distributed 
formly within the aggregate, one has Ni){d/2)^ ^ 

i.e., Rb ^ Nl^^d/2, and A7V x (47r/3)(6//2)3 ^ AttRIAR, 

implying AN ^ SA/"^ (AR/d). Finally, assuming that 
the thickness of the interfacial layer, separating the seg- 
regated cell regions, is AR = xd^ where 2 < x < 3, one 
obtains 



^2/3 



(18) 



Note that according to Eq. (18), as A/" ^ 00, i.e., for large 
aggregates, Smax approaches unity as N~'^^^ (assuming 
that Na and Njj are of the same order of magnitude). 

The time evolution of the sorting parameter, s = s{t), 
in our KMC simulation corresponding to case CI is shown 
in Fig. [?[ The insets represent snapshots of the sorting 
process taken at times indicated by the arrows. 

The sharp increase of s{t) at the beginning of the sim- 
ulation followed by a slow asymptotic approach to Smax 
indicates that there are at least two sorting time scales. 
Indeed, the entire time evolution of the sorting parameter 
can be well fitted with the double exponential 



sit) 



sie 



-t/ri 



S2e 



-t/T2 



(19) 



where Smax = 0-76 is in very good agreement with the 
theoretically estimated value 0.78 obtained from Eq. (18) 
for X = 2.5. The other fitting parameters in Eq. (19) are: 
ri = 1.4to, si = 0.27, rs = bSMo and S2 = 0.11. The 
shorter time scale ri corresponds to the local rearrange- 
ment (sorting) of cells leading to small clusters of same 
types of cells, while the longer time scale r2 describes 
the much slower engulfment process of the b cells by the 
a cells, a process that requires large displacements by a 
finite number of cells. 

Although the results of our KMC simulations appear 
to be in good qualitative agreement with experiments on 
cell sorting O [T3[ l28^ , a quantitative comparison, e.g., 
in terms of the time evolution of the sorting parameter, is 
not feasible because s{t) cannot be measured experimen- 
tally. Thus, there is no simple way to reliably calibrate 
the time unit to (which is related to the model parame- 
ter wq) used in the plot of s vs t/to in Fig. ItI However, 
s{t) can also be determined from CPD simulations, thus 
allowing for a quantitative comparisson between the two 
computer simulation methods. 



2. CPD simulations 

We have also used CPD simulations to investigate cell 
sorting corresponding to the three cases CI, C2 and C3 
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FIG. 7. Time dependence of the cell sorting parameter, 
s = s(t), corresponding to case CI described in the text, for 
both KMC (top) and CPD (bottom) simulations. The insets 
represent snapshots of half of the spherical aggregate taken 
at times indicated by arrows. 



described above. The initially spherical aggregate con- 
tained a random mixture of equal number Na = N^j = 
1, 000 of cells of type a and b. While the CPD parameters 



0.8 and 655 



Antra 



1.2 



were kept the same in all three simulations, the parame- 



ter Cab 



Anter 



had different values (similar to the ones 



used in the KMC simulations) for the three cases CI, 
C2 and C3 as listed in Table [ij The cell sorting patterns 
obtained at the end of the corresponding CPD simula- 
tions are shown in Fig. |6] As expected, these patterns 
are similar to the ones obtained in the KMC simulations. 

In order to quantify the degree of cell sorting in the 
CPD simulations by employing the cell sorting parameter 
5, defined through Eq. (17), we determined the position 
of a cell by the center of mass of the constituent CPs, 
and considered two cells to be neighbors if they were 
separated by a distance less than 3.25 a. For the CPD 
simulation corresponding to case CI, s{t) is shown Fig.JTJ 
Similarly to the KMC result, s{t) can be fitted well with 
the double exponential (19). Again, Smax = 0.68 is in 



good agreement with the theoretical prediction Eq. (18), 
i.e., 0.67 for x = 2.2 (or 0.63 for x = 2.5). The other 
fitting parameters in Eq. (19) are: ri = 0.68 to, r2 = 
103 to, si = 0.25 and S2 = 0.1. Note that while si 
and 52 have essentially the same values for both KMC 
and CPD simulations, the time constants ri and T2 are 
quite different, as the corresponding time units to are 
different in the two simulations. Moreover, the fact that, 
for similar model parameters, r2/ri = 41.8 in KMC is 
about twice as large as r2/ri = 21.7 in the corresponding 
CPD simulation indicates that the self diffusive motion 
of cells in KMC occurs much faster than in CPD. In 
other words, the multicellular system is more liquid-like 
in KMC than in CPD simulations. 



IV. CONCLUSIONS 

We have presented two (KMC and CPD) simulation 
methods and an analytic, continuum theoretical ap- 
proach to address structure formation by cell sorting and 
the fusion of contiguous multicellular spheroids. The the- 
oretical method was used to interpret the experimental 
results, and to test the viability and calibrate the model 
parameters of the KMC and CPD simulations. Our study 
was motivated by the need to quantify biomechanical 
properties of engineered tissue constructs, composed of 
compact tissues made of adhesive and motile cells and 
to predict their time evolution. The growing interest for 
understanding shape changes in such tissue constructs 
stems from their applications in tissue engineering in gen- 
eral and in the emergent field of 3D bioprinting in par- 
ticular [8 . 

The KMC method is based on a lattice representa- 
tion of the 3D tissue construct and dynamics is described 
in terms of rates associated with possible movements of 
cells. Similarly to previously employed MMC studies, the 
mixing pattern observed in KMC simulations disagrees 
with experiments. In both methods an elementary move 
consists in cells swapping positions with neighbors, which 
overestimates cell motility. Nevertheless, the time scale 
calibration in KMC makes the simulated time course re- 
alistic as far as the shape evolution of multicellular tissue 
constructs is concerned. 

The CPD method is based on modeling individual cells 
in a tissue construct as interacting CPs. The dynamics 
of the multicellular system are determined by integrat- 
ing the equations of motion for each CP. The CPD force 
field parameters are determined such that the time evolu- 
tion of the shape of the fusing spherical aggregates in the 
CPD simulation matches as closely as possible the exper- 
imental one (i.e., two touching spherical caps). Once the 
CPD model is calibrated, this can be used to simulate the 
shape evolution of arbitrary 3D multicellular constructs. 
It should be emphasized that in CPD (i.e., computer) 
units the calibrated CPD parameters (and therefore the 
outcome of a CPD simulation) are independent of the 
used cell type. However, the CPD units ([9| have spe- 
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cific values for different cell types. Thus, the CPD sim- 
ulations reported here can be applied as is to different 
cell types; the corresponding CPD time unit to should 
be determined in each case by equating Tgim ^ 540 to 
(see Sec. Ill A 2) with the experimental fusion time Texp- 
The reported CPD simulations provided a good de- 
scription for both fusion and cell sorting of multicellular 
spheroids. We found that CPD provides a more realistic 
description of complex multicellular structure formation 
than KMC. Indeed, the behavior of the studied multi- 
cellular systems in CPD simulations resembles to that 
of complex visco-elastic materials while in KMC simula- 
tions to that of viscous liquids. It is to be expected that 



by including more realistic features into the interaction of 
the CPs the accuracy of the CPD method can be further 
improved. 
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